Simulation of two dimensional quantum systems on an infinite lattice revisited: 
corner transfer matrix for tensor contraction 
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An extension of the projected entangled-pair states (PEPS) algorithm to infinite systems, known as 
the iPEPS algorithm, was recently proposed to compute the ground state of quantum systems on an 
infinite two-dimensional lattice. Here we investigate a modification of the iPEPS algorithm, where 
the environment is computed using the corner transfer matrix renormalization group (CTMRG) 
method, instead of using one-dimensional transfer matrix methods as in the original proposal. We 
describe a variant of the CTMRG that addresses different directions of the lattice independently, and 
use it combined with imaginary time evolution to compute the ground state of the two-dimensional 
quantum Ising model. Near criticality, the modified iPEPS algorithm is seen to provide a better 
estimation of the order parameter and correlators. 
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Understanding the emergent properties of many-body 
systems is one of the main goals of modern Physics. Pro- 
jected entangled-pair states (PEPS) [I] were recently pro- 
posed by Verstraete and Cirac to describe the ground 
state of finite, inhomogeneous quantum systems on a 2D 
lattice. A PEPS consists of a two-dimensional network of 
tensors whose coefficients are optimized so as to approx- 
imate the ground state of a local Hamiltonian. An ex- 
tension of the PEPS algorithm to infinite, homogeneous 
2D lattices, known as infinite PEPS (iPEPS) algorithm 
[2], was also subsequently developed. As a many-body 
ansatz, the infinite PEPS had already been discussed by 
Sierra and Martm-Delgado [o] under the name of vertex 
matrix product ansatz, and had been successfully used 
by Nishino and Okunishi [4], under the name of ten- 
sor product variational state, to evaluate the partition 
function of a 3D classical system. In the context of 2D 
quantum systems, a simplified version of the ansatz, with 
only three free parameters, had also been used in Ref. 
[5] prior to the iPEPS algorithm [2]. By considering an 
evolution in imaginary time, the iPEPS algorithm over- 
came the stability problems of previous proposals while 
allowing for the optimization of all the coefficients in 
the ansatz. Subsequently, ingenious approaches to op- 
timize homogeneous PEPS in finite systems with peri- 
odic boundary conditions have also been proposed, such 
as those in Refs. [(>] and [7], where the ansatz is often 
called tensor product state. So far, algorithms based on 
the PEPS formalism have provided accurate ground state 
properties of several models of spins and hard-core bosons 
[1, 2, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Importantly, they 
can address systems beyond the reach of quantum Monte 
Carlo, such as frustrated antiferromagnets [12, 11]. 

An approximation |^) to the ground state of a local 
Hamiltonian H with an infinite PEPS is typically ob- 
tained either by minimizing the expected value of the 
energy (^'jiJI^') or by simulating an evolution in imag- 
inary time l^*) « e~^'^|^'o), where |\E'o) is some initial 
state. In either case, the tensors that define the infinite 



PEPS are optimized iteratively and, in order to properly 
update a given tensor, one needs to compute its envi- 
ronment: a 2D tensor network that accounts for the rest 
of the tensors in the ansatz. Unfortunately, computing 
the environment is hard and an approximation scheme is 
required. In the case of a homogeneous system, several 
approximation schemes have been proposed: 

(i) ID transfer matrix approaches, such as the infinite 
time-evolving block decimation (iTEBD) approach [14] 
used in the original iPEPS algorithm [2]. 

(ii) 2D coarse- graining approaches based on the tensor 
entanglement renormalization group (TERG) [(i, 7]. 

(iii) Corner transfer matrix (CTM) approaches, such 
as the corner transfer matrix renormalization group 
(CTMRG) algorithm [J ]. 

In this work we explore a modification of the iPEPS 
algorithm where instead of computing the environment 
using the iTEBD approach, as originally proposed in [''], 
we use the CTMRG [ ]. 

The CTM formalism was originally derived by Bax- 
ter [ ] and later adapted by Nishino and Okunishi in 
the CTMRG [15] to numerically compute environments. 
The goal of this paper is to investigate the performance 
of CTM approaches to compute environments within the 
context of the iPEPS algorithm [ ] . Specifically, we first 
describe a versatile variant of the CTMRG, the direc- 
tional CTM approach, which addresses different direc- 
tions of the lattice separately, and use it in conjunction 
with imaginary time evolution to study the 2D quan- 
tum Ising model near criticality. Accurate estimates of 
the critical magnetic field and P exponent are obtained. 
Then we compare the results obtained using the original 
iPEPS algorithm (where the environment was computed 
using the iTEBD [14]) with this new version of the algo- 
rithm. The modified iPEPS algorithm is seen to converge 
signficantly faster to the ground state and provide a bet- 
ter characterization of the critical point [ ] . 

Let us consider an infinite PEPS for the state l^*) of 
an infinite 2D lattice C, which for concreteness we take 
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FIG. 1: (color online) Diagrammatic representation of (a) 
infinite PEPS tensor A with physical index s and bond in- 
dices u, r, d and I; (b) reduced tensor a; (c) infinite 2D tensor 
network £; (d) environment f for site r; (e) eight-tensor 
effective environment C/'^- 



FIG. 2: (color online) (a)-(d) Main steps of a left move: in- 
sertion, absorption and renormalizatlon; (e) the CTMs Ci, 
Ci and the half-row transfer matrix T4, are renormalized with 
isommetry Z; (f) eigenvalue decomposition for the sum of the 
squares of CTMs Ci and C4. 



to be a square lattice, with each site labeled by two in- 
tegers r — [x, y) and represented by a complex vector 
space V of dimension d. In the simplest scenario, the in- 
finite PEPS is characterized by a single tensor A that is 
repeated on all lattice sites. It has components As udir, 
where s labels a local (physical) basis of V (s = 1 , • • • , c? 
) and u, d, I, r are bond indices ranging from 1 to D, with 
D the bond dimension of the infinite PEPS, see Fig. 1(a). 
Let a denote the reduced tensor a = J2t=i^s ^ A*, 
with double bond indices such as u — {u,u') (see Fig. 
1(b)). Then the scalar product (^'|^') can be expressed 
as a two-dimensional network £ made of infinitely many 
copies of a, Fig. 1(c). [Notice that with a proper choice 
of tensor a, £ can also represent the partition function 
of a 2D classical statistical model] . The environment of 
site r, fl*^ = d£/da^^, is obtained from £ by removing 
the tensor a on site r, see Fig. 1(d). The goal of the 
CTMRG algorithm is to compute an approximation Q^^ 
to f ['T by finding the fixed point of four CTMs. This 
effective environment is given in terms of a small ten- 
sor network, G^^ = {Ci, Ti, C2, T2, C3, Tg, 6*4, r4}, where 
tensors Ci , C2 , C3 , C4 represent four CTMs (one for each 
corner), and tensors Ti,T2,r3,T4 represent two half- 
column and two half-row transfer matrices. Fig. 1(e). 

In the directional variant of the CTMRG used in this 
work, the eight tensors of G^^^ are updated according 
to four directional coarse- graining moves, namely left, 
right, up and down moves, which are iterated until the 
environment converges. Given an effective environment 
g[rl = {Cl, Ti, C2, T2, C3, T3, C4, T4}, a move, e.g. to the 
left, consists of the following three main steps. Fig. 2: 



(1) Insertion: insert a new column made of tensors Ti, 
a and T3 as in Fig. 2(b). 

(2) Absorption: contract tensors Ci and Ti, tensors C4 
and T^, and also tensors T4 and a, resulting in two new 
CTMs Cl and C4, and a new half-row transfer matrix 
fi, see Fig. 2(c). 

(3) Renormalization: Truncate the vertical indices of 

Cl, Tj^ and C4 by inserting the isometry Z, Z'^ Z — I. 
This produces renormalized CTM's C^ = Z^Ci, C^ = 
C4Z and half-row transfer matrix T4, Fig. 2(d)-(e). 

A proper choice of isometry Z in the renormalization 
step is of great importance. One possibility is to use 
the eigenvalue decomposition of the product of the four 
CTMs Cl, C2, C3, ^4 as in Ref. {V% Here we consider 
instead the eigenvalue decomposition of CiC| -I- C4C4 — 
ZDzZ'^ , Fig. 2(f), and use the isometry Z that results 
from keeping the entries of Z corresponding to the x 
largest eigenvalues of Dz- Similarly to Ref. [15], this 
isometry targets the CTMs of the effective environment 
instead of the wave function itself. The cost of imple- 
menting these steps scales with D and % as OiD^x'^. 
The net result is a new effective environment Q'^^ for 
site r given by tensors {C^ , Ti, C2, T2, C3, T3, C^, T^}, see 
Fig. 2(d). By composing the four moves of the direc- 
tional CTM we recover one iteration of CTMRG [lo]. 
The additional flexibility provided by individual moves 
can be used to accelerate convergence in a specific direc- 
tion, e.g. in highly anisotropic systems. In addition, the 
prescription used to compute the isometry Z is still valid 
-and produces stable results [ I T]- in the context of sim- 
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ulating imaginary time evolution described in this work. 
As with other similar methods [14, l^>, 19], an imme- 
diate application of the directional CTM is to compute 
expected values from 2D classical partition functions (re- 
sults not shown). 

In order to compute the ground state of 2D quantum 
models by simulating imaginary time evolution, we con- 
sider an infinite PEPS characterized by two tensors A and 
B. The optimization of tensors A and B proceeds in the 
same way as we proposed as part of the iPEPS algorithm 
[2] , but with the crucial difference that here the required 
environment for two contiguous sites is computed with 
the directional CTM instead of using ID transfer matrix 
(IDTM) techniques. The scalar product (^'1^') consists 
of an infinite 2D tensor network made of copies of the re- 
duced tensors a and b. We first consider the environment 
£;['-i,i-2,r3,r4] q£ ^ four-site unit cell, see Fig. 3(a), and ap- 
proximate it with an effective environment Cy[''i>''2.''3,''4] — 
{Ci, Tf,i, Tai, 6*2, T'a2, '7b2, ^3, Tfjs, Ta3, C4, Ta4, T^} made 
of twelve tensors, which are computed by iterating left, 
right, up and down moves. These directional moves are 
a natural adaptation to a four-tensor unit cell of those 
in Fig. 2. As shown in Fig. 3(b), the half-row and half- 
column transfer matrices Ti (i = 1,2,3,4) are replaced 
with pairs of half-row and half-colum transfer matrices 
Tai,Ti,i. This time, in order to implement e.g. a left 
move, two new columns are inserted in the system in 
step (1), see Fig. 3(c). For each of the inserted columns, 
we perform the absorption and renormalization steps (2) 
and (3). The renormalization step requires introducing 
an additional isometry W , which we compute in an anal- 
ogous way as isometry Z, see Fig. 3(e). As before, the 
cost of a move scales as 0{D^x^)- Finally, from a con- 
verged environment for the four-site unit cell, an effective 
environment for any pair of nearest neighbor sites is eas- 
ily obtained with an additional directional move. 

To demonstrate the performance of the approach, we 
have computed an infinite PEPS approximation to the 
ground state of the spin- 1/2 quantum Ising model on a 

transverse magnetic field, Hj{\) = —J2{rr'}'^^^'^^^^ ~ 
A^-crlf^, by simulating an imaginary time evolution. 
The simulation proceeds as in Ref. [:], but we use the 
directional CTM to obtain the effective environment at 
each step of the imaginary time evolution. This evolution 
is performed with decreasing time steps St ranging from 
10"""^ to 10~^, and until convergence of local observables 
and two-point correlators is attained. 

Fig. 4 shows the order parameter = (^'jo-jl^') as 
a function of the transverse magnetic field A, for {D, x) 
equal to (2,20) and (3,30), where the value of x is cho- 
sen so that the results are converged with respect to this 
parameter. Remarkably, an infinite PEPS with bond di- 
mension D = 3 already produces results within less than 
a percent from the best quantum Monte Carlo estimates 
for the critical magnetic field A*^'-^ w 3.044 and criti- 




FIG. 3: (color online) (a) Environment of the four-site unit 
cell; (b) twelve-tensor effective environment; (c) two new 
columns are inserted, and absorbed towards the left and 
renornalized individually. The diagram shows the contrac- 
tion leading to Ci,C4, Ta4 and Tt4 when absorbing the first 
column, and also to the CTMs Qi and Qa; (d) two isometries 
Z and W are used to obtain the renormalized half-row trans- 
fer matrices T^4 and T54; (e) eigenvalue decomposition for the 
sum of squares of CTMs Qi and Q4,. 





iPEPS with 


iPEPS with 


TERG [ ] 




directional CTM 


ITEBD [2] 




Ac D=2 


3.08 


3.10 


3.08 


D=3 


3.04 


3.06 




13 D=2 


0.333 


0.346 


0.333 


D=3 


0.328 


0.332 





TABLE I: Critical point Ac and exponent (3 for the 2D quan- 
tum Ising model as estimated by the new and old versions of 
the iPEPS algorithm, as well as the TERG (for a finite lattice 
of up to 2^ X 2^* spins). For reference, the quantum Monte 
Carlo estimation is \^'^ « 3.044 and p'^''^ ^ 0.327 [:>!)]. 



cal exponent for the order parameter p-^'^ 0.327 [20], 
namely with relative errors ~ 0.1% and « 0.3% respec- 
tively. Table I contains a comparison with results ob- 
tained with the original version of the iPEPS algorithm 
and with the TERG algorithm for large systems [ ]. 

It is particularly instructive to compare the perfor- 
mance of the original and present versions of the iPEPS 
algorithm, since they are both based on imaginary time 
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FIG. 4: (color online) Order parameter as a function of 
the transverse field A, computed with the directional CTM 
approach. Lines are a guide to the eye. The lower-left inset 
shows a log plot (in natural logarithms) of rriz versus |A — Ac|, 
including our estimates for Ac and /3. The continuous lines 
show the linear fits. The upper-right inset shows a compari- 
son close to criticality with the results from Ref.[ ] using the 
original version of the iPEPS algorithm, which used iTEBD 
(dashed lines). Results correspond to {D,x) equal to (2,20) 
and (3,30). 

evolution and only differ in how the two-site environment 
is computed: by means of the iTEBD and directional 
CTM approaches, respectively. One finds that when com- 
puting environments using the directional CTM, a sig- 
nificantly better infinite PEPS approximation to near- 
critical ground states is obtained, leading to a more ac- 
curate characterization of the quantum phase transition. 
As shown in Fig. 5, the resulting infinite PEPS also dis- 
plays stronger correlators S'zz(0 = (^'lo'i'^cri^^''^"'' j^*) — 

However, further comparison of results involving also 
other spin models reveals that, away from the quantum 
critical point, both the directional CTM and the iTEBD 
approaches yield equivalent accuracies for ground state 
properties. In particular, both versions of the iPEPS 
algorithm are equally suited to study first order phase 
transitions, a task for which they are particularly suc- 
cessful [9, 13]. Indistinguishable results are also obtained 
in models with long or infinite correlation lengths, such as 
the Heisenberg antiferromagnet on a square lattice [1 ], 
when using a bond dimension D that is too small to of- 
fer a proper approximation to the ground state (in Ref. 
[13], a bond dimension D = 5 still produces a sponta- 
neous magnetization that is off by 10%). It seems, there- 
fore, that computing environments with a CTM approach 
leads to better results when the following two require- 
ments are simultaneously met: (i) the ground state must 
have a long correlation length (e.g. near a quantum crit- 
ical point), and (ii) the bond dimension D must be suf- 
ficiently large that the ansatz can in principle properly 
approximate the ground state. 

The correlators in Fig. 5 also show that numerical cal- 
culations with infinite PEPS introduce an artificial finite 




FIG. 5: (color online) Log plot of the correlator SzziJ) (in 
base 10), as computed with the original and present versions 
of the iPEPS algorithm, namely using the iTEBD (dashed 
lines) and the directional CTM (solid lines) approaches. Lines 
are a guide to the eye. Our results are for the (D, x) pairs 
(2, 20) (for A = 3.09) and (3, 30) (for A = 3.04). 

correlation length at criticality. This is a consequence of 
the truncation in parameter x in the calculations of the 
effective environments. This truncation effect is similar 
to the one discussed in Ref. [21] and could in principle 
be analyzed in a systematic way. 

To summarize, PEPS are a valuable ansatz to approx- 
imate the ground state of 2D lattice models, as previ- 
ously demonstrated by several authors [1, 2, 5, G, 7, 8, 
9, 10, 11, 12, 13] with a variety of systems of interact- 
ing spins and hard-core bosons, including frustrated spins 
[12, 13] that cannot be addressed with quantum Monte 
Carlo techniques. Several methods have been proposed 
to optimize the PEPS. The main factor limiting the ac- 
curacy of the results is the bond dimension D. However, 
for a fixed value of the bond dimension D, the quality of 
the approximation may also depend on the method used 
to optimize the ansatz. In this work we have investi- 
gated a modification of the iPEPS algorithm, where the 
iTEBD of the original proposal has been replaced with 
a directional CTM, a variant of CTMRG [ ]. The new 
version of the algorithm provides a significantly better 
description of the ground state near criticality. 
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